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The behaviour of elastic structures undergoing large deformations is the result of the competition 
between confining conditions, self- avoidance and elasticity. This combination of multiple phenomena 
creates a geometrical frustration that leads to complex fold patterns. By studying the case of a 
rod confined isotropically into a disk, we show that the emergence of the complexity is associated 
with a well defined underlying statistical measure that determines the energy distribution of sub- 
elements, "branches", of the rod. This result suggests that branches act as the "microscopic" degrees 
of freedom laying the foundations for a statistical mechanical theory of this athermal and amorphous 
system. 

^ ! PACS numbers: 46.32,+x, 46.65.+g, 64.70.qd 

<5-> ■ I. INTRODUCTION 

m 

Classical equilibrium statistical mechanics stands as a major cornerstone of modern physics. Tools issued from 
this theory have been instrumental in rationalizing a huge number of seemingly unrelated physical situations ranging 
from phase transitions in atomistic systems to the behaviour of polymers. This is possible because the size of the 
fundamental components of those systems is sufficiently small so that thermal fluctuations allow the degrees of freedom 
to span all the possible configurations through an ergodic exploration of the energy landscape. For macroscopic systems 
thermal agitation becomes negligible and, while those systems may be mechanically stable, they fall out of equilibrium 
with respect to thermodynamics. A usual example where grav itational energy completely dominates thermal effects 
concerns the physics and mechanics of granular assemblies [32j . Because of this mismatch of energy scales and of the 
f^j \ presence of dissipative interactions between the grains, it had long been believed that those systems could never be 
reconciliated with the basic principles underlying classical statistical physics. However, it has been recognized in the 
last decade that granular materials do in fact share many properties with thermal, although very slowly evolving, 
systems such as molecular glasses [ll|, [24|. Indeed, analogies between the slow relaxation of granular materials and 
the glassy phenomenology of supercooled liquids has led to the now famous "Jamming diagram" proposed by Liu and 
Nagel [41] . While this analogy is not free from criticism, it does open up the possibility that macroscopic systems could 
still be described using a combination of the tools of statistical physics and some new out-of-equilibrium concepts. 
In fact it has recently been proposed that the rheology of athermal amorphous systems is characterized by well 
defined statistical distributions that are amenable to a probabilistic treatment similar to that of classical statistical 
' mechanics @, 0, 0, [H, |H EH . An important new concept often encountered in those descriptions is the presence of 
an "effective" temperature that is different from the bath temperature and that depends on the driving mechanism 
or on the details of the coupling between the thermostat and the system [3, HH EH EH • 

Here we are interested in another kind of athermal macroscopic systems: elastic structures. Due to the geometrical 
frustration created by the interaction between confinement, self- avoidance and elasticity, it is expected that confined 
elastic structures should display a complex behaviour. Indeed crumpling a sheet of paper is a typical everyday illus- 
tration of the phenomena that we want to study. A rapid visual inspection of the numerous valleys and mountains 
present after unfolding a piece of crumpled paper often leaves one with the impression of fascinating albeit extremely 
complex fold patterns [57j . Accordingly, the folding phenomena are associated with a rich class of crumpling phenom- 
ena which belong to a wider class of interfacial deformation phenomena. Both deterministic and random folding of 
thin materials are of noteworthy importance to many branches of science and industry. Examples range from DNA 
packaging inside virus capsids [Ho| and polymerized membranes at the microscale [l6|, [5l| to folded engineering mate- 
rials and geological formations at the macroscale, including insect wings [17] or leaves in buds [2l|, l34j . They usually 
consist of thin sheets or rods constrained to undergo large deformations. Because of their biological and technological 
importance, the properties of randomly folded thin materials are now the subject of increasingly growing attention. 
Because self-avoidance and nonlinear deformations make the description of fully crumpled or folded materials very 
difficult, earlier studies focused on the identification of the elementary generic features displayed by an elastic surface 
that has to accommodate a geometrical mismatch reducing its accessible volume. It was found that those individual 
structures consist of sharp vertices (or developable cones) [8|, Il8l420l [36| and linear ridges [26|, 1421444 [58| and their 
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properties are now rather well established. These singularities are conceptually similar to the dislocations, cracks, 
necks or shear zones that are created when a continuum media is forced at the large scales and localizes the stress 
in those tiny regions that dissipate the energy at the smaller length scales. In general, these singularities can either 
be moving or be quenched depending on their type and on the mechanical properties of the material [HI, H3] . Their 
interactions are usually carried through a long ranged elastic field either instantaneously or with retardation via wave 
propagation. In the case of the packing of flexible structures, almost nothing is known about how these interactions 
lead to the complex fold patterns that are observed at the surface of a piece of crumpled paper. This situation is 
similar to the one encountered in statistical physics where, given the knowledge of the microscopic interactions, one 
tries to bridge the gap and rationalize the macroscopic behaviours [l5|, [33| . 

The question that motivated the study presented here is: "Are there any general statistical properties associated 
with the crumpling and folding phenomenology that is observed when an elastic structure is confined in an environment 
smaller than its rest size?" So far there has been some degree of incompatibility between crumpling experiments and 
numerical simulations. For one thing, most simulated sheets have been fully elastic [35|, [37], l54i l56| whereas the ones 
used in experiments have been made of elasto-plastic materials such as Mylar, paper, aluminium foil and even layers 
of cream jl|, 0-0, OH [3l|, l38l440[ lifSl |49| . It is not clear to what extent these analyses could make the difficult 
distinction between elastic and elasto-plastic behaviours [55]. On the other hand, the packing of elastic objects clearly 
depends on the dimensionality of the object (d) and of the container (D > d). The configurational properties are 
totally different if the dimension of the container satisfies D > d + 2 from the case D = d+ 1 due to the rigid constraint 
of the inter-penetrability. 

Here we provide a case study in order to answer our original ques tion. We are studying the statistical properties 
of an elastic rod (d = 1) that is confined in a disk (D = 2) [~L3ll25l l27l |28[ l37l l53| . In this geometry, there are no 
singularities and the most important constraint is of geometrical origin. Using a comparison between two widely 
different implementations of rod compaction we suggest that there exists an underlying statistical measure that 
describes the folding process. After describing the systems that we studied (sheet pulling experiment in section ITTT1 and 
minimal numerical simulations in section HV|) and pointing out their different confining and geometrical conditions, we 
show a compared statistical analysis of several quantities. Similarities between stochastic observables can be detected 
by comparing their probability density function (pdf) and this is the method that we employed here. The particular 
pdf 's that we observed during this study are presented and put in a more general context in section [TTJ Our main result 
is that while the topological (section [V)) and the geometrical (section [VT]) properties are indeed different, the energy of 
some sub-parts of the rod (referred to as "branches" afterwards) display identical statistical properties (section IYlI|) . 
The energy of the branches turns out to be distributed according to a Gamma law that reduces to a Boltzmann 
distribution for high energies. In section IVIII) we argue that branches correspond to the relevant "microscopic" 
degrees of freedom that define the statistical mechanics of folding. Finally we propose that the spiral configurations 
that define the ground state of the system may be viewed as a condensed phase with all the branches lying on top of 
each other. 



II. STATISTICAL DISTRIBUTIONS 



Since we will be extensively using probability density functions (pdf) in the remaining part of the article, it is worth 
reminding the reader of the expressions of the pdf's that will come up later. In order to be general, we denote by x 
the random variable in this section but that will of course be replaced by some physical observables such as curvature 
or energy when we come to the description of our results. In addition to defining our notation for the parameters of 
the pdf's, this short list also allows us to summarize results from other related studies. 

• Exponential distribution: 

Pe( x ) = ;; ex P (-;;); (!) 

where \i is the mean of x. Such pdf's have been found before in measurements of local 3d curvatures [lOj and 
number of intersecting ridges [l[ in unfolded sheets of paper. They also appeared to describe the lengths of folds 
in a 3d crumpled sheet [56j . Exponential distributions are usually associated with the presence of uncorrelated 
events that are distributed randomly. 

• Log-normal distribution: 
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FIG. 1: a) Schematic representation of the sheet pulling experiment, showing the radius of the sheet r, its thickness h, the radius 
of the hole R and the control parameter Z . b) Typical configurations observed after a cross- sectional cut right below the hole. 
Note that these very different shapes were obtained by following the same protocol (compaction rate e — 19). c) Corresponding 
configurations after the image analysis procedure (see subsection \IIIB}) . Black circles represent the Y-shaped junction points 
where multi-branched stacks converge (merge or separate). Branches in self-contact within a stack make the thickness of these 
stacks proportional to the number of individual branches they carry. 

where \i and a 2 are respectively the mean and the variance of \nx. Such pdf's have been found to describe the 
length of plastic linelike ridges in experiments p], [Toj| as well as in numerical simulations of crumpled sheets [54], 
|56| . Log- normal distributions are usually associated with random events that are occurring hierarchically. 

• Gamma distribution: 

padHx) = m^ exp \x) ; (3) 

where x is the mean of x and a is the exponent of the Gamma distribution. Such pdf's are known to describe 
the lengths of bent segments of highly confined sheets [54[ . Gamma distributions are usually associated with the 
presence of correlations in otherwise randomly distributed events. When a < 1, there is a power- law divergence 
for small values of x. Notice that for large values of x, a Gamma distribution reduces to a simple exponential 
distribution. 

We will keep the same notation for all the parameters throughout the rest of the article. Also, we will not re- write 
the expressions for the pdf's and refer the reader to this section whenever one of these pdf's is encountered in the 
following sections. The goal of our study is to identify what is the relevant observable x leading to a universal pdf in 
the context of confined elastic rods. The robustness of the pdf's is tested by comparing between experimental work 
(set-up described in section HIT]) and "model" numerical simulations (algorithm described in section |IV|) . 

III. SHEET PULLING EXPERIMENTS 
A. Experimental set-up 

The experiment has been already described in [l3|, [25|. A circular polyester sheet of radius r ~ 30 cm and 
thickness h ~ 0.1 mm is pulled by its center through a smaller circular rigid hole of radius R ~ 2 cm, as illustrated in 
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L(m) 


h/L 


r/L 


R/L 


Zm/ L 


£ 


C 


B(J) 


1 


1.9 


2.6 10" 5 


8.9 10~ 3 


0.18 


0.16 


19 


15% 


710" 5 


2 


1.9 


2.6 10~ 5 


1.2 10~ 2 


0.17 


0.16 


14 


8% 


710~ 5 


3 


1.2 


2.1 10" 4 


2.2 10" 2 


0.18 


0.16 


8 


9% 


110" 3 


4 


1.0 


5.0 10~ 5 


1.6 10" 2 


0.33 


0.16 


10 


4% 


710" 5 



TABLE I: Experimental parameters: dimensions, compaction rates and material properties of the polyester sheets. 




FIG. 2: Zoom on a typical configuration detailing the successive steps of image processing, from left to right. Raw image. 
Pixels of blue values larger than a given threshold but of red and green values smaller than another given threshold are kept. 
Thresholded image. Empty spaces of area larger than Smin are kept. Skeleton of the binary image and junction points. 

Fig. [TJ The Young's modulus of the sheets has been measured as E = 5 GPa, its density is 1.4 g/cm 3 . The bending 
rigidity has been measured as B = 7 10 -5 J for sets of experiments 1,2,4 and B = 1 10 -3 J for set 3. In this specific 
set-up, the sheet undergoes preferentially bending deformation rather than stretching, in order to minimize its elastic 
energy. Thus a self-affine conical shape is expected and observed, so that the whole 3d shape of the sheet is prescribed 
by the shape of one cross-section. Any cross-section at a distance z from the cone tip, draws a virtual rod of thickness 
ft and length L ~ 2ttz, compacted in a circle of surface S = tt(Rz/Z) 2 , Z being the pulling distance between the cone 
tip and the plane hole (Fig. [1]). In this article, we are specifically interested in the compaction of such a virtual rod, 
rather than the compaction of the whole sheet. Lengths, surfaces and volumes measured in a cross-section obviously 
depend on its position z, but become independent on z, when their units are non-dimensionalized by using the total 
length of the rod L. In practice, the cross-section is observed in the hole plane z = Z, and only once during the 
pulling experiment for Z = Z m ~ 30 cm (see Table [Tj). For each experiment, a rod of non-dimensionalized length 1 is 
compacted in a circle of non-dimensionalized surface S/L 2 : 

S/L 2 = R 2 /4/Zl=e 2 /A, (4) 

defining the compaction rate £, the ratio of rod length over hole perimeter: 

e = L/(27rR). (5) 

A second parameter characterizing the compaction rate, but in terms of higher dimensionality, is defined as the ratio 
of the sheet volume over its conical envelope : 

C = 3he 2 /r. (6) 

Note that the sheet dimensions (ft, r) appear in the expression of the 3d-compaction rate C, but not of the 2d- 
compaction rate e. In this article, the experiments are typically characterized by compaction rates e of the order of 
10 and C of 10%. The values of experimental parameters L, ft, r, R, Z m , £, C, E and B are reported in Table HI 

The compaction mechanism in the present experiment is geometrically controlled by imposing the total available 
surface S for a rod of length L. The use of a rigid hole constraint can be viewed as a hard- wall repulsion potential 
V(r) acting on the rod such that V(r) = for r < R and V(r) = oo elsewhere. For fixed control parameters (E, ft, 
e and C), several pulling experiments have been performed in order to obtain various configurations corresponding 
to different energy minima and to analyse them in a statistical way. We realized four sets of experiments at constant 
control parameters (Table UJ). Some typical examples of folded rod obtained are shown in Fig. [TJ 



B. Image analysis procedure 

Because of technical difficulties [25|, we resorted to a hot wire cutting tool to obtain cross-sections in the hole 
plane. With great care, one obtains neat cuts without perturbing the configuration, and inks them in blue to reinforce 
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contrast of the cut edge with back receding surfaces. The cross-section is digitized with a scanner at a resolution of 
50 pixels/mm, which yields 5 pixels for the thickness of the thinner sheets (Fig. [2j). A thresholding based on RGB 
values results in a binary image: only pixels of blue values larger than a given threshold but of red and green values 
smaller than another given threshold are kept, allowing to distinguish the cut edge from back receding surfaces and 
to remove noise from the raw image (Fig. [2]). Then, empty spaces of surface area larger than s m i n = (10h) 2 are kept, 
allowing to remove light noise from the binary image (Fig. [2j). Finally, the binary image is skeletonized -reduced 
to a one pixel thick skeleton, without redundant kinks- (Fig. [2]). Junction points are then defined as pixels with at 
least 3 neighbours (Fig. [1] and Fig. [2]). Two neighbouring junction points delimit a stack of branches in close contact, 
separated by distances smaller than the image resolution. 

The next step is to determine the number of branches n^ r in each of the M stacks. The conservation of the number 
of branches at each junction point yields 2M/3 equations, because 3 stacks intersect at each junction point. This 
requires to distinguish converging stacks in 2 opposite packs: merging and separating stacks. To this aim, the vectors 
tangent to the stacks locally to the junction point are determined. Then, the straight line perpendicular to one of 
these vectors is taken as a reference axis, that delimits two semi-planes: stacks located in the same (resp. opposite) 
semi-plane are in the same (resp. opposite) pack. The result is systematically checked to be independent on the choice 
of the reference axis. The remaining M/3 equations are found from the thickness of the stacks in the binary image as 
follows. The heating by the cutting tool thickens a stack nonlinear ly, which was calibrated by separately cutting stacks 
of sheets. We keep the M/3 stacks with the best estimation of the thickness as given by the calibration (in general, 
thin stacks). The solution of the linear M x M system yields the number of branches in each stack. We reopened a 
few configurations (5 per set of experiments) and checked by counting the number of branches in each stack: we found 
no error for sets 2, 3 and 4 (corresponding to e < 14, C < 9%), and an error of ±1 for a part of the thicker stacks 
(20% of the stacks) in the more compact set 1 (corresponding to e = 19, C = 15%). These errors are small thanks 
to the fact that the number of branches is an integer, and to the successive steps of image processing. The central 
position of stacks is thus deduced from the skeleton and smoothed out by slide- averages along the stack trajectory. 
The local curvature is measured through a parabolic fit on size 5 h, after coordinates translation and rotation in 
the tangent frame. Image processing as described above allow to detect branches of length t > tmin = 1 ram and 
voids of surface s > s m i n = (why. 



IV. NUMERICAL SIMULATIONS 



A. Energy functional 



We consider an elastic rod of bending rigidity B and total length L. Its configurations are represented by a 2d 
vector R(s) parametrized by the arclength s G [0,L]. Contrary to the experiment described in section Hill the global 
constraint of compaction is introduced by plunging the rod into an external quadratic potential. In that case, instead of 
a hard- wall geometry, the external field acts as a body force attracting the rod towards the minimum of the potential 
located at the origin R = 0. Note that imposing a stronger constraint such as R a with exponent a ^> 2 (which 
would be closer to the experimental compaction potential) makes the minimization algorithm significantly slower. We 
will see later that this is not a limitation as we do not intend to precisely mimic the experiment in the numerical 
simulations but rather to extract some robust observations common to both systems. The strength of this confining 
potential can be varied through a control parameter A. Once embedded in this confining force-field, the shape of the 
rod is prescribed by the competition between two effects: While its bending rigidity tends to keep it straight, the 
rod responds to the external confining force by buckling and developing folds. Since we will consider infinitely thin 
rods in the simulations, the resulting deformations can only be of pure bending and the rod is unstretchable: Its total 
length L remains a conserved quantity. Elasticity theory shows that the energetical cost for bending deformations is 
proportional to the square of the rod curvature. In addition to the bending energy, there is another source of energy 
associated with the confinement. The total energy of the confined rod can be written as: 

B f L /d 2 R\ 2 A f L 
E = — / ( ) ds H — / R 2 ds + "Hard-core repulsion". (7) 



2 J V ds 2 J 2 J0 

The physical constraint of self-avoidance is implemented through a discontinuous hard-core interaction. In the fol- 
lowing, the spatial variables R and s are rescaled (now denoted R and s) by the rod length L and the total energy 
E by B/2L. This yields the total dimensionless energy E: 

E = J Q ( ~Jf ] d5 + A J ^ 2<i5 + "Hard-core repulsion" , (8) 
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FIG. 3: Simulations. Typical random initial configurations used as input in the numerical simulations. These are generated 
by imposing a white noise on their curvature according to Eq. [TUjj) , with noise intensity D = 750. 

with the dimensionless control parameter A = XL 4 /B. In order to compute this energy numerically, the rod is 
discretized into TV segments of constant length L/N. Derivatives for the bending energy are determined via finite- 
differences and integrals are computed with usual trapezoidal rules. Our goal is to explore the energetical landscape by 
minimizing Eq. (|5J). An obvious constraint that has to be satisfied during the minimization is self- avoidance (physical 
objects cannot cross themselves). This constraint is encoded in the "hard-core repulsion" term of Eq. (jSJ). Whenever 
a configuration contains at least one self-intersection, its energy is set to infinity (a very large number in practice). 
Otherwise only the first two terms (bending and confinement) contribute to the total energy for configurations free 
of any self-intersections. 

Avoiding self-intersections in the energy minimization is a delicate operation because it is a non-local interaction. 
Regions of the rod that are far away in the rest state (straight rod) become very close when the available area decreases 
and folds start to appear. Since the location and the nature (localized or extended) of these contact regions cannot 
be known beforehand, the detection and treatment of self-contact areas is numerically expensive. Because the rod is 
discretized into a connected polyline of TV segments, looking for self-intersections is a procedure that usually involves 
testing each pair of segment and therefore grows as N 2 . However here, we can take advantage of the particular 
(connected) geometry of the problem and lighten this procedure. Instead of testing each pair of segments, the idea is 
to restrict our search to a limited set of segments that are more likely to contain self-intersections. This is possible by 
designing a variation of the brute force N 2 method that keeps track of the distance between segments. Once a pair of 
segment has been tested we also determine their distance in units of segment length L/N. This information is then 
used to determine how many of the following segments can be ignored. These segments are skipped because they are 
too far to generate an intersection with the tested segment (even in the worst case scenario where they would come 
straight back towards it). Therefore the actual number of tested segments depends on the input curve. We found this 
input sensitive algorithm to run much faster than the simple N 2 method and as fast as more elaborate techniques 
(sweep-line for example) in computational geometry, at least for the moderate values of N ~ 300 used here. 

B. Minimization procedure 

As mentioned above, the self- avoidance property introduces discontinuities in parameter space by building up 
infinite barriers between different attraction basins. In general, the question of finding extrema of discontinuous 
functions is very difficult, because one cannot use familiar procedures such as gradient methods. We found Powell's 
algorithm to be convenient in our minimization problem. It is a derivative-free procedure whose search directions in 
parameter space are updated after each iteration, finally generating non-interfering (or conjugate) search directions. 
The numerical protocol is now described step by step: 

1. Initial configurations. A catalog of random initial conditions is constructed by the introduction of a white 
noise of amplitude D on the curvature of the rod such that : 

<R») = 0, (9) 
<R"(«i)R"(a 2 )> = 2DS( Sl -s 2 ), (10) 

where the average is taken over each realization of the rod and 5 is the Dirac distribution function. Once a value 
for D has been set, we can iteratively generate multiple initial shapes. The configurations' center of gravity is 
translated back to the origin. Some typical configurations used as random initial shapes are shown in Fig. [3l 
The largest accessible value of D is the one which does not generate self-intersections. 
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FIG. 4: Simulations. Typical local minima obtained with A = 7 10 5 and D = 750. The corresponding mean compaction ratio 
is (e) = 2.3. 

2. Introduction of the confinement. The parameter A controls the strength of the confining potential. Going 
instantaneously from A = for the construction of initial configurations, to a non-zero value to take into 
account the confinement, can be seen as a quenching mechanism. This is because the confinement constraint is 
instantaneously turned on to its desired value without taking intermediate steps. All of the results presented 
here were obtained by following this procedure. We also tested an annealing process by increasing A very slowly, 
about which we will say few words in the conclusion. 

3. Search for local minima. In practice we ran Powell's algorithm 10 times for each initial configuration. 
We verify at each iteration that configurations do not contain any self-intersections. If they do, we set the 
energy of such configurations to a very large number, so that these configurations are immediately rejected 
by Powell's algorithm. The set of search directions is regularly re-initialized to a set of random directions in 
order to optimize our span of local minima. Running the minimization procedure many times also improves 
numerical convergence to these local minima. One should note that while the hard core repulsion procedure 
ensures self- avoidance, it raises problems when parts of the rod are aligned. Indeed our numerical procedure 
does not allow sliding of self-contact areas. Therefore we (temporarily) modeled self-contacts as a nematic 
interaction of the form -ELi f _ CO ntact = u sin 2 a where u is a dimensionless parameter and a is the angle between 
touching segments [l5|, [33]. Typically, we used u in the range 50-150 but its precise value did not affect the 
resulting configurations. This geometric self-interaction has a physical interpretation as a self-excluded volume 
resembling what was introduced by Onsager in the context of polymers. In our case, it destabilises tightly 
clamped configurations by allowing branches in self-contact to self-align and slide along one another. This 
nematic term is included only in 2 early algorithm runs (out of 10) and ignored otherwise (and in particular in 
the last run). 

After going through this procedure a number of times with different initial conditions but under the same compaction 
conditions (noise intensity D and control parameter A), a wide variety of different possible configurations is found, 
some typical examples of which are shown in Fig. [U The results presented here are with A = 7 10 5 and D = 750 for 
a total of over 250 different realizations. 



C. Simulations vs experiments 

As previously noticed, the global compaction constraint is of different nature in simulations and experiments: 
compaction is mechanically controlled in simulations through A, the intensity of the quadratic potential exerted on 
the rod, whereas it is geometrically controlled in experiments, through the available size R for the rod. The pressure 
exerted on the quasi- Id rod in experiments is related to the force necessary to pull the 2d sheet through the hole. 
This raises the question of the parameter common to experiments and simulations, relevant for the description of the 
compaction strength. On one hand, pressure is not trivially accessible in experiments; on the other hand, the size of 
the occupied surface can be easily characterized for each configuration in simulations by its radius of gyration R g : 




allowing to compute the compaction rate e, written in Eq. (|5j), by replacing R by R g . The radius of gyration R g 
has been computed for all the numerical configurations. Its probability density, presented in Fig. [5j features a 
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FIG. 5: Simulations. Distribution of the radius of gyration R g non-dimensionalized by the rod length L, of the 252 configura- 
tions with A = 7 10 5 and D = 750. 



FIG. 6: Definition of elementary branches. Parts of the rod delimited by two junction points (•) define the elementary 
"branches" . Some of them are in close contact and are delimited by the same junction points, defining a multi-branched 
stack. Experiments: configuration from set 3 at compaction rate e — 8, with a snapshot showing a multi-branched stack. 
Simulations: the dashed circle represents the mean radius of gyration (R g ) at e — 2.3. 

sharp peak at (R s ) ~ 0.05L, showing that imposing a confining potential and noise amplitude in the simulations, 
results to indirectly imposing the size of the surface occupied by the rod. In these early stages of folding that are 
accessible numerically, we define the geometric compaction rate as e = L/ (27ri? g ), and thus have on average for all 
configurations (e) = 2.3. 

To summarize, the values of e studied in this paper are collected for experiments and simulations respectively. 
Four sets of experiments are realized at fixed compaction rate e = 8, 10, 14 and 19 (see Table |T]). Even if the 
compaction rate is not a priori controlled in the simulations, it is indirectly imposed to the mean value (e) = 2.3. 
Despite the different nature of compaction, geometrically and mechanically controlled in experiments and simulations 
respectively, implying a different expression of the rod energy (as it has been suggested previously and it will be 
shown later), and despite the different values of compaction rates e ~ 10 in experiments and e ~ 2 in simulations 
(visible when comparing Fig. [T] and Fig. 2]), both systems share the property that these are governed by elasticity and 
self- avoidance. The main issue of this study, allowed by the investigation of these two systems, is to point out the 
general characteristics of folding and to identify the ones that are system dependent. Two types of averages will be 
used in this study: x refers to the mean of the variable x over one configuration, whereas (x) refers to the ensemble 
average of x over the whole set of configurations realized at constant control parameters. Due to better statistical 
averaging, most of the results presented here were obtained by ensemble statistics. 

A crude observation of examples of folded rods shows that, under constant control parameters, a wide variety of 
configurations are obtained in both the experiments (Fig. [1]) as well as the numerical simulations (Fig. H]). This is an 
interesting first remark confirming that we do indeed span a large volume of phase space. Although the shape of a 
folded rod looks very complicated, elementary parts of the rod, delimited by two neighbouring junction points (where 
the rod is locally in self-contact) can be identified. These elementary supports are referred to as "branches" , and could 
be relevant candidates for a "microscopic" definition and parametrisation of the macroscopic folded configuration. 
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(N br ) 


(nbr) 


(s)/L 2 


(V~s/P) 


(Z) / L 
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380 
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6.6 10~ 6 


0.19 


2.5 10" 3 
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80 
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40 
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TABLE II: Averages over all realizations of the number of multi-branched stacks Ni, number of branches Nb r , number of 
branches per stack nbr = Nb r /Ni, surface area s, shape ratio y/s/p of elementary tiling voids and length of branches i. 
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FIG. 7: Experiments, a.) Mean number of branches per stack Hbr as function of the total number of stacks Ni. Continuous 
(resp. dashed) lines correspond to Eq. \13\) (resp. Eq. \15\) ). h) Square root of the mean area of voids as a function of 
the mean length of branches I (both quantities are non-dimensionalized by L). The continuous line corresponds to a linear fit 
of slope 1.1. c) Illustration of the exponential distributions of the number of branches per stack nb r for the different sets of 
experiments. 

From a mechanical point of view, branches are the natural elements, on which elastic equations are defined given 
the boundary conditions at their extremities. Fig. [6] shows two folded configurations coming from experimental and 
numerical sets, where junction points are shown by circles, that delimit branches. Branches that are in close contact 
and delimited by the same junction points define a multi-branched stack. We start our data analysis in section IVl 
and section [VTl by looking at the topological and geometrical properties. While their behaviour are interesting, it is 
in fact the energetical properties studied in section IVTI1 that turn out to be the desired general variable. 

V. GENERAL PROPERTIES OF THE NETWORK 
A. Stacking 

Here, we are interested in the number of different elementary units in a folded configuration: junction points, 
voids, branches, multi-branched stacks, and to their relative dependencies. Junction points, multi-branched stacks 
and number of branches in close contact define respectively the nodes, the links and the values attributed to each link 
of a network Therefore, we start by investigating the topological properties associated with the network formed 
in the confined rods. 

Experiments. A folded configuration is characterized by a number Ni of multi-branched stacks. Each stack is 
formed by a number n^r of multiple branches, the sum of which is A^ r over the configuration. A natural question is 
how these network properties Ni and vary together? The values of (N{) and (A^ r ), reported in Table HH increase 
roughly with the compaction rate. Fig. [T^i shows the mean number of branches per stack = N^/Ni, as a function 
of the total number of stacks TV/, for all experiments of all sets. Inside a single set, plotted points seem to follow a 
systematic trend (Fig. [7^), showing that A^^ r and Ni are correlated. 

Exactly 3 multi-branched (or links) intersect at 99% of the junction points, the value 3 being the minimal possible by 
construction. This allows to relate the total numbers of multi-branched stacks Ni and of voids N v in a configuration: 
Ni « 3N n /2 w 3N V . For a closed rod, the number of junction points (or nodes of the folded network) N n is equal 
to the number of voids 2N V (or cells of the network). By definition, N v ~ S/s ~ Ni/3 and 7V^ r ~ L/I, where s and 
I are respectively the mean surface area of voids and the mean length of branches in a configuration. Indeed the 



10 



total length of the rod L is divided in the N^ r branches and the total available surface S is tiled by the N v voids. 
Moreover, L and S are related through the control parameter e = Lj This allows to write the following 

relation between Nb r and Ni: 

N br = 2 ~§e^VNi , (12) 

through the parameter £ and an a priori configuration dependent factor y/jj/l. It is equivalent for the mean number 
of branches inside one stack averaged over configuration to: 

Nbr 1 
nb r = — — = — -^e— -= . (13) 

Formula ([12]) and (fT3|) are valid for a single experiment, are these still verified for a whole set of experiments? Can one 
write a relation between A^ r and Ni through quantities averaged over sets of experiments (instead of configurations), 
or through scalar factors? This raises the issue of the variations of y/^/l inside a set of experiments: is this ratio 
varying with Ni and control parameters el To this aim, we plot the square root of the mean voids area as a 
function of the mean branches length I for all experiments in Fig. 03. Even if these values change by factors 2 
inside a set and by a factor 5 between different sets, these appear to be proportional through the constant value 1.1, 
independent of the precise configuration and of the compaction parameters, as shown by the linear fit in Fig. 03. 
So one can write from Eqs. ([T2]) and (fT3j) : Nb r oc e\/Wi and n& r oc e/y/Ni, that are plotted in continuous lines in 
Fig. [7K, super-imposed on experimental data. We will see in the following sub-section that another relation between 
Nbr and Ni can be predicted. There is still an issue with the understanding of the universal value of \/I/7, whatever 
the configuration, the compaction rate and the size of the sheet. 

We found that the number of branches inside one stack rib r follows an exponential distribution as shown on Fig. 0c. 
For the different sets of experiments, the mean (n& r ) of the exponential distribution varies from 3 to 5 for increasing 
e (Tab. [TTJ) . Inside the sets of data, the relative fluctuations {ribr / (nbr) — 1) are only of the order of a few percents 
meaning that (nbr) — (nbr}- A precise analysis of the pdf p(jibr) would require much more statistics. Nevertheless, it 
is tempting to make the assumption that n& r ~ (n& r ). This implies that n& r is not correlated with other configuration 
dependent variables such as Ni. Therefore, this reasoning leads to a prediction of the total number of branches Nb r 
and the average number of branches per stack nbr per configuration, as functions of the number of stacks Ni given by: 

N br = (n hr ) N L , (14) 
n br = (nbr) • (15) 

These linear predictions are plotted as dahes lines in Fig. [T^i for a comparison with Eqs. (fT2j) and (fT3|) and with the 
experimental data. It turns out that both predictions ([T2]) - ([T3]) and (jHJ-JTSJ) describe well the data. This suggests 
that Ni and Nb r are in fact selected by the validity of both assumptions. The intersection of the two relations 
corresponds to: 

m vm ^ 2 _ 2 ^ (16) 

e 

At the level of accuracy accessible in our experiments, the agreement is not so bad even though there seems to be a 
small dependence with e. 

Simulations. As the rod is open, the surface it encloses is not well-defined and not exactly constant from realiza- 
tions to realizations. The analysis presented for the experimental results cannot hold for the simulations. However, 
Fig. [8k shows the number of branches N hr as a function of the number of multi-branched stacks Ni for all the numerical 
simulations: these verify the relation N hr ~ 2N^ the distribution of N hr /N is indeed strongly peaked around 2 as 
shown on Fig. [Hfr. Several configurations are characterized by the same pair (A^ br , Nf). 

Because the compaction rate achieved in the numerical simulations, (s) = 2.3, is much smaller than those of the 
experiments, the number of branches in self-contact is only rarely more than 3. Generically links contain just 2 
branches. This explains the previous macroscopic relation N hT c± 2N U shown in Fig. [U 



B. Tiling 

We focus here on the properties of the elements of tiling of the total available surface, such as their perimeter, 
surface and shape. Note that these elements correspond not only to loops (defined as an elementary "closed" part 
of the rod in [28|), but also to the voids that are delimited by such touching loops. Does the spatial tiling of the 
available surface contain any information on how the folding might have happenned? 
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FIG. 8: Simulations, a Number of branches Nb r as a function of the number of multi-branched stacks Ni in the simulations. 
Continuous line corresponds to a line of slope 2. b Distribution of the ratio Nb r /Ni. 
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FIG. 9: Experiments, el) Distribution of the area s of the voids. Continuous lines correspond to fits by log-normal distributions, 
h) Distribution of the aspect ratio p(y/s/p) where p is the perimeter of the voids. It is well described by a Gaussian distribution. 
The values 0.28 for a circular shape and 0.22 for a pinched closed elastic shape are drawn as vertical lines, c) Distribution of 
the number of links m surrounding a void. 



Experiments. We found that the shapes as well as the sizes of the voids show a wide variability over 2 or 3 
orders of magnitude. The distribution of the area s of the voids shown in Fig. [9^ is well described by log-normal 
distributions. The mean value (s), reported in Table HH decreases with the compaction rate. The shape of the voids 
can be characterized by the ratio of the square root of its surface over its perimeter y/s/p, as done in [28j. This aspect 
ratio is shown in Fig. [9]b. The distributions are well described by an unique gaussian distribution independent on the 
control parameters: 

*"<■)- vE"^) • <17) 

with p = 0.19 and a 2 = 14 10~ 4 . One can also look at the number of junction points belonging to each void 
contour, that is also the number of edges in terms of multi-branched stacks surrounding the void. One should note 
that a junction point that belongs to a void contour does not correspond necessarily to a vertex where geometric 
properties (local tangent or local curvature) of the two adjacent stacks change discontinuously but rather to a change 
of the number of branches in the stacks. The voids do not have a constant number of edges n\ as shown by the 
distribution in Fig. [9fc. The paucity of our data does not allow us to discriminate between Gamma, Log-normal or 
exponential distributions as they all describe well the data n\ > 4. 

Simulations. As already discussed, the simulations are not realized at constant available surface. We defined the 
surface S occupied by the open rod as the sum of surfaces of elementary voids; this distribution is shown in Fig. fTUh. 
The pdf p(S) is peaked around an average value of (S) = 4 10 _3 L 2 that has the same magnitude as that expected 
for a perfectly circular shape tt (R g ) 2 = 8 10 _3 L 2 . On the other hand, the distribution p(s) of the surface of the 
elementary voids changes over 3 orders of magnitude and is closer to a power law distribution as shown on Fig. [TUb. 
The shape of this distribution is unchanged regardless of whether the surfaces of the voids are normalized by the 
occupied surface S or by the mean void surface s of the configuration. Fig. [TUb shows the pdf of the aspect ratio 
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FIG. 10: Simulations, a.) Distribution of the surface area S occupied by the rod. h) Distribution of the area s of the elemtary 
voids, c) Distribution of the aspect ratio p(y/s/p) where p is the perimeter of the voids. It is well described by a Gaussian 
distribution. The values 0.28 for a circular shape and 0.22 for a pinched closed elastic shape are drawn as vertical lines. These 
statistics are based on 1,435 elementary voids. 

y/s/p of voids: it is widely distributed of mean equal to 0.17, and well described by the same gaussian distribution 
than for the experiments. 

C. Summary 

To summarize, we observed the difference of values of Ni and 7V& r between experiments (~ 50 and ~ 200) respectively 
and simulations (~ 10 and ~ 20) due to the different compaction rates e (~ 10 and ~ 2). The number of elementary 
units growing with the compaction rate underlies the increasing complexity of the configurations. Several predictions 
(Eqs. (fT2]) - (fT5]) ) are obtained on the relative variations of Ni and either from geometrical relations or from 
observations at the scale of multi-branched stacks. They both appear to be validated by the experimental data. 
Eqs. (fT2]) - (fT3|) result of the independence of the ratio y/^/l on the configuration, the compaction rate and the size of 
the sheet (Fig. [Tfr): there is still an issue with the understanding of this constant value. Note the large range and 
values accessible by (Fig. [7]c), the number of branches in close contact forming a multi-branched stack: it would 
be interesting to study what controls the value of n& r inside one multi-branched stack. However, whatever the degree 
of complexity of the folded configurations, some statistical characteristics stay unchanged, namely the shape ratio 
and the number of edges of the elementary voids. The aspect ratio of the voids also appears to be independent on 
the compaction parameters e (varying from 2 to 20 between simulations and experiments) and on the compaction 
potential (hard- wall or quadratic). The average value of the aspect ratio of the voids is 0.18. A perfectly circular 
shape provides an upper bound of y/s/p = 0.28. For comparison, this ratio is y/s/p = 0.22 for a clamped elastic 
rod as can easily be shown by solving the Euler's elastica equations. The values we report here can be explained 
by considering that the loops are confined: their perimeter will not change but the area will be reduced due to the 
compression. The experiments reported in [28| show values y/s/p ~ 0.17-0.20, very close to 0.22, consistent with their 
injection method that tends to create a layered loops geometry. Finally, the distributions of voids surface have been 
found to be Log- normal in the experiments. This hints that voids are generated through a hierarchical process. A 
succession of bifurcations such as the one predicted in [l3| leads to folding events that successively "break" voids. A 
new fragmentation event concerns all the previously existing fragments. However, this pdf is closer to a power law 
distribution in the case of the numerical simulations. Due to the small range of accessible compaction ratio, it is 
difficult to assess the relevance of this observation. 



VI. GEOMETRICAL PROPERTIES OF BRANCHES 



In the following, we will characterize the geometry of the folded rod at two scales. On the one hand, we will 
investigate the geometry of branches, namely theirs lengths and mean curvatures. On the other hand, we will look at 
the local curvature along the rod. However, let us first define two sub-systems in the experimental set-up. 
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TABLE III: "Microscopic" properties averaged over experiments from the same set. Geometry: mean curvature K m and length 
of branches £. Energy: Total energy E, energy of branches e, and parameters of the Gamma distributions a e and Xe (resp. as 
and xe) for branches energy (resp. total energy), according to Eq. (0). 




FIG. 11: Experiments. Definition of two sub-systems, a.) Branches that have one of their extremities in direct contact 
with the container wall are colored in red to separate them from the other branches that form the bulk of the system in green, 
b) Super-imposition of the two sub-systems for experiments of set 1. Note that because there is no well defined boundaries in 
the numerical simulations, this subdivision is unique to the experimental system. 



A. Two sub-systems in the experiments. 

The branches composing the folded rod are differentiated into two sub-systems: branches with (resp. without) an 
extremity in contact with the container, which we will refer to as periphery (resp. bulk). Fig. [TT1 shows branches in 
the bulk and in the periphery on one example of configuration and for all super-imposed configurations from set 1. 
Statistical properties of branches can be studied either in the whole system, or separately between the periphery and 
bulk sub-systems. We found that there is systematically more branches located at the periphery than in the bulk. 
On average, the periphery (resp. the bulk) is composed of 60% (resp. 40%) of branches. Defining branches in the 
periphery as the ones for which both extremities are in contact with the wall would have resulted in a slightly different 
outcome. However, our choice of definition allows to prevent any boundary effect on the bulk properties. For example, 
we found that the average number of branches per stack is the same in the bulk and periphery; thus it is uniform in 
the whole system. 

B. Length of branches 

Experiments. The distribution of the length of branches is found to be exponentially distributed in both sub- 
systems as shown in Fig. [T2h-b. It appears that the mean length is slightly larger for branches at the periphery than 
in the bulk. The ratio of mean lengths between the two sub-systems (£ pe riph) / (hulk) is between 1 and 1.6 depending 
on the set of experiments. Using the fact that approximately 40% (resp. 60%) of the branches are located in the bulk 
(resp. at the periphery), the mean length in the whole system verifies {£) ~ 0.4 {hulk) + 0.6 (iperiph)- The resulting 
pdf in the whole system is therefore exponential as shown in Fig. [T2b . Quantitatively, the average length of branches 
(£), reported in Table ITTT1 decreases with the compaction rate. The values for branches in the sub-systems are written 
in the legend of Fig. [T2l 

Note that the existence of a minimal value £min, related to experimental detection limits, may introduce a small 
shift of the average (£) towards larger values, in comparison with the characteristic decreasing exponential length fi: 
(£) > \i. Strictly, experimental biased data are distributed according to: (£) exp(— £/ji)/ii 2 . However, as the value 
£min is very small in comparison with (£), this later well approximates \i. This method has the advantage to be 
insensitive to the choice of bins for the construction of the pdf, contrary to a simple minded fit. 
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FIG. 12: Experiments. Geometrical properties. Distribution of the length £ of branches a): In the bulk sub-system, and 
exponential distributions of means \i — 2.1 10 _3 ; 2.710 -3 , 7.3 10 -3 and 4.3 10 -3 for sets 1-4. b): In the periphery sub-system, 
with means \i — 2.9 10 -3 , 4.0 10 _3 ; 12 10 -3 and 4.5 10 -3 . c): For all branches in the whole system, with means \i — 2.5 10 _3 ; 
3.410 _3 ; 10 10 -3 and 4.4 10 -3 . Vertical lines show the minimal length £min = 1 mm experimentally detected. 
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FIG. 13: Simulations. Geometrical properties. Probability distribution p(£) of the length of branches £ in log-lin and log-log 
scales. The continuous line represents a Gamma-law with parameters: a = 0.4 and % = 0.12. Control parameters are: A = 7 10 5 ; 
D = 750 and there are 4676 branches. 



Simulations. The distribution of the length of the branches follows a Gamma law (Eq. (|3j)) with shape parameter 
a smaller than 1: a = 0.4 and x = 0-12 as shown in Fig. [T3j This observation is important because it indicates an 
accumulation of branches of small lengths, in comparison with the case of a pure exponential pdf which was observed 
in the experiments. 



C. Curvature of branches 



Several definitions are possible for the characterization of the branch curvature: 

k± = / K(s)ds/£ , = / \n(s)\ds/£ , and n m — \k±\ , (18) 
Jo Jo 

A 2 

with k(s) = ^j-irS the local curvature along the branch at curvilinear abscissa s. When a branch does not exhibit an 
inflexion point the three definitions above are equivalent, with n± having possibly a different sign. The statistical 
study of mean curvature of branches ft m allows for a representative sampling of all of the values of local curvature 

Experiments. We found that the whole rod has an increasing number of inflexion points as e increases. However, 
because inflexion points ususally coincide with junction points connecting adjacent branches, the local curvature does 
not change change sign at the branch level. Indeed, only less than 1% of the branches have an inflexion point 
where k+ ^ n m . Thus, the statistical distributions of k± appear to be exactly symmetric around 0, so that the 
pdf of \k±\ is equivalent to that of n m . The pdf of the mean curvatures in the bulk are approximately distributed 
according to exponential laws as shown in Fig. fT4h. By contrast, the same pdf in the periphery is characterized by 
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FIG. 14: Experiments. Geometrical properties. Probability distributions p(ft m ) of mean curvature of branches K m for sets 
1-4- a: In the bulk sub-system and exponential distributions of mean (him) = 304, 230, 93, 160. b: In the periphery sub-system: 
(nm) — 350, 290, 110, 230. c: For all branches in the whole system: (ft m ) = 330, 260, 100, 200. The exponential pdf in h-cdo 
not have the same mean than the data, but comes from a. 
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FIG. 15: Simulations. Geometrical properties, a: Probability distribution of the absolute curvature per branch ft+ and 
Gamma distribution. Parameters for the Gamma-law are: a = 2.5 and x — 14.8. b: Probability distribution of the curvature 
of branches k± with a normal distribution (dashed line) with: p — and a — 39, and an exponential distribution (continuous 
line) with p — 30. Statistical test: \ 2 — 93 for the exponential distribution and 160 for the normal law. Data is averaged over 
4,676 branches. 



a symmetric and sharp peak around the curvature of the hole Kh which is clearly visible in the inset of Fig. [T4b. 
The presence of this peak means that, even though both pdf's are well described by exponential laws, their means 
are significantly different. As a result, the pdf p(ft m ) of mean curvatures ft m for all branches in the whole system, 
verifying p ~ OAp^ik + 0.6p per i P h, is also qualitatively a global exponential law plus a symmetric peak around the 
hole curvature (Fig. [Tib). 

In addition, we found that lengths £ and mean curvature radius ft" 1 are slightly correlated (linear correlation 
coefficient r c± 0.45). Note that describing the branches as circles arcs, of perimeter smaller than for the corresponding 
whole circle, gives the inequality t < 27r/ft m , that is experimentally verified. 

Simulations. The distribution of ft + is well described by a Gamma distribution with parameters a = 2.5 and x — 
14.8, as shown in Fig. [T5a . The average over all configurations is (ft + ) = a\ = 37/L. As the linear scales plot 
shows (and corroborated by the shape parameter a greater than 1), there is a power law drop-off for small values. 
Otherwise, away from the peak the probability falls off with an exponential behaviour. This peak is interpreted as 
the result of "effective" walls located at distance R g from the center which leads to an accumulation of curvature 
at =37- 1/R g = 20. 

As the compaction ratio is much smaller than in the experiments, many branches have inflexion points in the 
numerical simulations. Therefore n+ and k± have very different statistical distributions. Note that whereas ft + gives 
direct information on the curved state of the considered branch, k,± alone is not obvious to interpret. A branch can 
have a small value of ft±, because it is either straight or because it contains an inflexion point. Therefore, a branch 
of curvature fc + is bounded by \k±\ < Fig. fT5b shows the pdf of k±. While it is symmetric around 0, it can 
equally be fitted by an exponential distibution of mean p = 30, or by a gaussian distribution of variance a = 39 (and 
mean 0). 
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D. Summary 

We interpret the emergence of exponential distributions for the branches length £ (pure exponential in the exper- 
iments and just an exponential tail in the numerics) as the random division of the whole rod. The length averages 
can be related to the number of branches. For each configuration we have £ = L/N^. One expects that this is still 
approximately valid for averages over all experiments of a same set, leading to {£) /L = (Nb r )~ (non-dimensionalized 
length). Thus the larger number of branches A^ r observed systematically for larger compaction rates £, explains 
the different quantitative exponential distributions. Another indication of the stochastic formation of the branches 
comes from the well representation of curvature statistics by the sampling made on branches. One common feature 
to experiments and simulations on branches curvatures statistics is their exponential distributions -pure exponential 
in bulk for experiments and for n± in simulations, or exponential tail for k+ in simulations-. These exponential 
distributions of 2d curvatures can be compared with the exponential distributions found in [lOj , for 3d curvatures in 
unfolded crumpled sheets of paper. 

After having looked at the geometrical properties, we will now turn our attention to the energetical characteristics 
of both experimental and numerical systems. 



VII. ENERGETICAL PROPERTIES 



A. Expression of energy 

In experiments, because of the self-similar conical shape of the folded sheet, related to pure bending deformations, 
its elastic energy is given by: 

BZ , f r \ f L /d 2 R N 2 



^-"■uj/. w) is - (19) 

with on the one hand, the bending energy B, a logarithmic prefactor based on the ratio of the sheet radius r and a 
cut-off length R c ~ 10 mm, and on the other hand, R and L the position of the rod and its length, observed in the 
cross-section at distance Z from the cone tip. By rescaling E s by the characteristic length scale due to the conical 
shape, Zln (r/i? c ), we obtain the energy of the rod per unit of transversal length: 

B f L /d 2 R N 2 



This expression for elastic energy per unit of length is the same than for elastic energy of the numerical rod (Eq. [7]), 
except that B in experiments has different units than B in simulations (energy and energy times length respectively) . 
We should mention that Eq. (fT9]) is exact only in the regime of pure elastic deformations. It turns out that for 
the highest compaction rates studied here, some configurations sustained a few localized plastic deformations. This 
happens when the absolute value of the local curvature \k\ goes beyond the plastic threshold. Through independent 
experiments, we measured it as n c = 0.54 mm -1 for sets 1, 2 and 4 and n c =0.24 mm -1 for set 3. In order to account 
for the subsequent plastic softening of the sheets in these areas, we substituted the quadratic dependance of E s on 
the curvature by a more appropriate linear dependance of the form k c (2k — k c ). This substitution does not affect the 
distributions, but only the total energy through extreme events. 

In simulations and contrary to the experiments, an important distinction has to be made right away about the 
expression of the energy of the rod in the numerical simulations. There are 2 independent terms contributing to the 
total energy: the bending energy and the confinement energy due to the quadratic potential. Because the numerical 
simulations are based on an energy minimizing principle, the total energy of the rod has already been discussed in 
Section HV( and is rewritten here (omitting the hard-core repulsion term which is for "equilibrium" configurations): 

B S L /d 2 R\ 2 . A ' L 



E =iJA^) ds+ 2L K2d '' (21) 



In non-dimensionalized form, the total energy of the rod can thus be rewritten as : 



experiments/numerics numerics only 
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FIG. 16: Experiments. Probability distributions of the energy of the branches, Gamma and Bose-Einstein distributions for 
sets 1-4. a: In the sub-system bulk: (e) = 290, 260, 95, 170, a e = 0.17, 0.23, 0.41, 0.24 and \e = 1700, 1100, 230, 720. b: In 
the sub-system periphery: (e) = 370, 407, 160, 260, a e = 0.17, 0.20, 0.44, 0.37 and Xe = 2200, 2100, 350, 710. c: In the whole 
system: (e) = 330, 340, 130, 220, a e = 0.16, 0.19, 0.41, 0.31 and Xe = 2000, 1800, 330, 720. Lower figures: Same as Fig.\M 
but represented in log-lin scale to better see the large exponential fall off of the distribution. 

where the first term (bending) is common to both the experiments and the numerical simulations, while the second one 
(confinement) is only present in the numerical simulations. In the following, we will be interested in two different levels 
of description (branch and whole rod scales) and Vt represents the domain of interest: "microscopic" (Q = branch) 
or macroscopic (ft = whole rod). Notice that in order to distinguish between the energy of branches and that of the 
whole rod, the notation e will refer to branches energy and E to the whole rod energy. 

B. Microscopic level (Q — branch) 

Experiments. For all the four sets of experiments, the distribution of the energy of branches p(e) is well described 
by Gamma laws whose parameters are given in the legend of Fig. [T6J As the local curvature of a branch is smoothly 
distributed around k m /k+ (section |VlJ), we verified that the energy of a branch is roughly equal to the estimation ik^ 
(the linear correlation coefficient is « 0.6). Taking into account the rough linear correlation between £ and /c" 1 , it 
implies a roughly linear correlation between k m and e. 

Simulations. Because the total energy of branches is made up of 2 independent contributions, it is interesting 
to first analyze them separately. It turns out however that both contributions, bending and confinement energies, 
are distributed according to a Gamma distribution with similar parameters (Fig. ITTa-b). This is a confirmation that 
the numerical procedure indeed converges to true minima of the energy functional. Summing up the bending and 
confinement energy yields the distribution of the energy of branches e (Fig. [TTb ) which is also well described by a 
Gamma law distribution. 

C. Macroscopic description (Q = whole rod) 

Instead of looking at the level of branches (microscopic), let us now focus on the total energy distributions of the 
whole rod (macroscopic). 

Experiments. The pdf of the non-dimensionalized total energy E of the rod seems asymmetric. In view of 
the comparison with the distribution of the branch energy, we also compare it to Gamma distributions as shown 
on Fig. HU Besides, Fig. fTDh shows that, individually, there is a significant correlation between the total energy of 
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FIG. 17: Simulations, a Probability distribution of the bending energy of branches. The continuous line represents a Gamma- 
law with parameters: a = 0.35 and x = 340. b Probability distribution of the confinenement energy of branches. The continuous 
line represents a Gamma-law with parameters: a — 0.23 and x = 492. c Probability distribution of the total energy of branches. 
The continuous line represents a Gamma-law with parameters: a = 0.31 and x — 728. Other parameters are: A = 710 5 ; 
D — 750 and there is 4, 676 branches. 
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FIG. 18: Distributions of total energy of the rod and comparison with Gamma distributions whose parameters (a, x) are - a-b) 
Experiments. (10.8, 1.2 10 4 ) ; (12.7,300), (15.9,40), (5.72, 1.3 10 4 ) and c) Simulations. (a = 61,x = 70). 

a folded configuration and its number of links/branches (linear correlation coefficient ~ 0.75), which confirms the 
central role of branches/links. 

Simulations. Just like we discussed for the microscopic description, the total energy is the sum of the bending 
energy and the confinement energy. Again, the distributions are slightly asymmetric. They can be described by 
Gamma distributions with similar parameters: bending {a = 27, \ = 81} an d confinement {a = 21, x = 100}. The 
sum of these two distributions yields again a Gamma distribution {a = 61, x = 70} for the total energy of the rod as 
shown on Fig. [T8h . On the other hand, Fig. [T9b shows that for the small compaction rates obtained numerically there 
is no clear correlation between the energy of the rod and the number of branches that it contains. We do however 
verify that the mean of the total energy is given on average by the average number of branches times their average 
energy: (E) = (N hr ) (e). 

VIII. CONCLUSION 

Despite some important differences beween the two systems (nature of confinement, compaction rates achieved, 
geometrical properties...), the major outcome of this comparison is the rather universal shape of the probability 
distribution of the energy of individual branches. In both cases, this probability density is well described by a 
divergent Gamma distribution: it has an exponential decay at large energies, and turns into a power law divergence 
for small energies (a < 1). This exponential fall off of the distribution is reminiscent of Boltzmann's law making it 
possible to extract a characteristic energy scale. The parameter x i n the Gamma distributions can then be considered 
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FIG. 19: To£a/ energy versus topological quantities, a) Experiments. Total energy E of the rod as a function of the total 
number of branches Nb r and the total number of stacks Ni for all experiments, b) Simulations. Total energy of the rod 
(normalized by (E) ) vs. the number of links (normalized by (Nu nks ) ). 
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TABLE IV: Test of the branches as the independent variables. While the experimental results suggest a convergence to a possible 
thermodynamic limit, the paucity of the numerical data cannot be interpreted in a simple way, see discussion in the text. 



as "effective" temperature characterizing the system. While we will come back below to the implications of this 
result, we can already say that it supports the notion of branches as the basic components of the folded rods (just 
like particles are in a real gas), interacting together through contacts forces at boundaries. At the other end of the 
spectrum, we do notice an accumulation of branches with very small energy (soft power law divergence) hinting at 
some underlying energy condensation process. 

The sum of N "uncorrelated" random variables distributed with a Gamma distribution of parameters {a^,x} is 
another Gamma distribution with parameters {a = Na^x\. Therefore the ratio a/a^ is a measure of the number 
of degrees of freedom present in the system. Assuming branches are the basic elements composing the whole rod, 
we would expect this ratio to be related to the number of multi-branched stacks. Table HVl shows that the expected 
relation beween a^/a e and the number of links is only crudely observed for the compaction rates achieved here. The 
experimental results show however that the higher the confinement and the better the agreement. This trend indicates 
that due to the small number of branches (particles) present in the configurations there are still non-negligible finite 
size effects introducing non-trivial correlations. As the number of branches increase, the system gets closer to the 
thermodynamic limit and the discrepancy indeed shrinks down a little. As the shape parameter a of the Gamma 
distributions, the probability density becomes more and more gaussian as a consequence of the central limit theorem. 
On the other hand, the numerical results completely fail this analysis. One possible interpretation is that, contrary to 
the experiments, one would expect branches rather than links to be the individual degrees of freedom in the numerical 
simulations. This is because branches are determined to be in self-contact whenever their distance is smaller than 
some small cut-off distance. This gives thick stack of branches more freedom to wiggle around and therefore behave 
somewhat more independently than in the experiments where there are true self-contact areas. In that case, one 
should perhaps think in terms of a sum of variables with Ni that are independent and N^ r — Ni that are exactly 
correlated. Our data can only hint in that direction but does not allow us any definite conclusion regarding this issue. 
In any case, the simple minded analysis we presented here is still quite interesting and would deserve more work to 
understand if one really finally converges to the thermodynamic limit as the number of branches tend to infinity. 

Bose-Einstein condensation and annealing. Another well-known statistical distribution that has the same 
asymptotic behaviours as Gamma laws is the celebrated Bose-Einstein distribution. Using the analogy between 
branches and gas-particles, we compared our data with this distibution at null chemical potential because the number 
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FIG. 20: Experiments. Averaged values of the energy of one branch e as a function of the number of branches in the 
multi-branched stack n& r . 
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FIG. 21: Simulations. Spiral configurations as the condensed phase of the system where all the branches are on top of each 
other. 



of branches (particles) is not pre-determined in either numerical simulations nor experiments: 

*> - 55^FT- < 23 > 

The agreement is qualitatively as good as the one obtained by using Gamma distributions. For high confinements 
(in experiments), it is easy to see from the experimental pictures that many branches tend to accumulate close 
to each other, giving rise to thick stacks of heavily populated branches. This phenomenon can be attributed to 
a transition from disordered configurations (isotropic phase) to ordered configurations (nematic phase) which was 
predicted theoretically. In addition, Fig. [20] reveals that these highly degenerate stacks (multiple branches close to 
of each other) do indeed carry very little elastic energy. It is interesting to note that this behaviour holds for low 
confinements as well. Because of their lower compaction rates, the configurations obtained numerically only rarely 
present stacks made up of more than 2 branches. However the Bose-Einstein distribution describes equally well the 
data. Since branches are not subjected to any exclusion principle (except for their vanishingly small thickness that 
prevents self-intersection), it is possible for them to be in the exact same state and effectively behave as integer 
spin particles explaining the observed good agreement of our data with the Bose-Einstein distribution. While this 
comparison may seem surprising, we will give below a few more arguments justifying our position. 

We saw that in addition to their exponential tail, the probability distributions of the branches display a power law 
divergence for small energies. By closely looking at the experimental images one can see that most of the branches 
tend to accumulate close to each other creating thick stacks made up of several branches. While this may not seem so 
obvious in the numerical simulations, one has to remember that they were obtained by quenching the rod to relatively 
much smaller values of compaction. Indeed if instead of the quenching mechanism proposed here, we increase very 
slowly the control parameter A (annealing), a different behavior is observed. We are no longer able to explore a wide 
phase space but instead converge to the true ground state. A typical shape that is obtained by this process is the 
spiral pattern presented in Fig. [2TJ We expect this phenomenon of branch condensation to the exact same state to 
hold and even increase for higher confinements. This is quite reminiscent of what happens during the condensation of 
integer spin particles (Bose-Einstein condensation) where particles are free to accumulate on the same energy levels. 
The analog of the chemical potential in our situation would be the number of branches present in the rod. 
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In summary, we have described two different systems (experiment and numerical simulation) to study the statistical 
properties of confined rods. The geometrical properties of the two systems are different due to the difference in the way 
confinement was applied. We do find however that the statistical properties of the bending energy of the branches 
follow the same distribution. Furthermore the parameters of this Gamma distribution share the same properties 
indicating a common behavior. 

Coming back to the question posed in the introduction, we can now answer to the existence of a statistical measure 
in this system. The energy of the branches is in fact the relevant internal variable. Moreover its distribution can 
be approximated over a wide range of energies by a Boltzmann distribution weighted by a characteristic energy 
reminiscent of the concept of "effective" temperature in the context of granular rheology. This system could be used 
to test some recent theorems related to non-equilibrium statistical physics in a context different from granular (and 
colloid) matter. An advantage of this system is that it is possible to measure the bending energy of branches which 
could be analogous to particles in a real gas. 
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